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We study dynamical and energetic instabilities in the transport properties of Bloch waves for 
atomic multi-component Bose-Einstein condensates in optical lattices in the tight-binding limit. We 
obtain stability criteria analytically, as a function of superfluid velocities and interaction parameters, 
in several cases for two-component and spinor condensates. In the two-species case we find that 
the presence of the other condensate component can stabilize the superfluid flow of an otherwise 
unstable condensate and that the free space dynamical miscibility condition of the two species can 
be reversed by tuning the superfluid flow velocities. In spin-1 condensates, we find the steady-state 
Bloch wave solutions and characterize their stability criteria. We find generally more regions of 
dynamical instability arise for the polar than for the ferromagnetic solutions. In the presence of 
magnetic Zeeman shifts, we find a richer variety of condensate solutions and find that the linear 
Zeeman shift can stabilize the superfluid flow in several cases of interest. 
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I. INTRODUCTION 

There has been considerable recent interest in the dy- 
namical properties of atomic Bose-Einstein condensates 
(BECs) in _optical lattice potentials [l|, 0, @, i, i, 0, i, 
i, llOl, nil, [il [H, [ii, 111 . In optical lattices the non- 
linear mean-field interaction of the BEC may give rise 
to dynamical and energetic instabilities in the transport 
properties of the atoms. For a single-component conden- 
sate, when the center-of-mass (CM) velocity reaches a 
critical value, the BEC dynamics become unstable result- 
ing in an abrupt stop of the transport of the atom cloud 
in the lattice. Such a superfluid to insulator transition 
has a classical nature and it can be described using the 
Gross-Pitaevskii (CP) mean-field models [13, E H H • 
In the dynamically unstable regime small initial pertur- 
bations around a moving solution grow exponentially in 
time, resulting in the randomization of the relative phases 
between atoms in adjacent lattice sites. The dynami- 
cal transition to inhibited atom transport was experi- 
mentally observed in the classical regime [2, 5, 6, ItI, [i| 
and experimental methods to characterize both the dy- 
namical and energetic instabilities of moving condensates 
have been developed Inhibition of transport was 

also observed in the presence of large quantum fluctu- 
ations using strongly confined narrow atom tubes ^] . In 
a confined system with enhanced quantum fluctuations 
the sharp classical transition is smeared out [13, El, , 
resulting in a gradually increasing friction in the atom 
transport. Due to the broadening of the velocity distribu- 
tion of the atoms, even at low velocities a non- negligible 
atom population occupies the dynamically unstable high 
velocity region of the corresponding classical system, gen- 
erating in the shallow lattice limit the friction ^u] . 

Despite this work on single-component condensates, 
there have been relatively few studies of dynamics of 
multi-component BECs in optical lattices. Due to nearly 



equal trapping potentials of different Zeeman sub-levels 
(for example F = l,mi? = —1 and F = 2, = 1 in 
^^Na and ^^Rb), it is possible to create long-lived two 
component BECs, forming an effective spin- 1/2 system. 
These have especially long lifetimes in ^^Rb due to a for- 
tuitous cancelation of scattering lengths [21]. This ad- 
ditional degree of freedom has been utilized to study an 
interesting array of effects in both Bose-condensed and 
non-condensed cold Bose systems, including phase sepa- 
ration_[22[ , optically-induced shock waves [23] , spin waves 
[H, llgToverlapping "^^K-^^Rb BEC mixturesJHK smn 
squeezing [27], and vector soliton structures [281, [29l, |3Q|| . 
Experimental work on two-component BECs in optical 
lattices, from the viewpoint of quantum logic gates, was 
reported in [31]. There has also been a recent experi- 
mental realization a two-species ^^K-^^Rb Bose mixture 
in an optical lattice [32]. 

Alternatively, in dipole traps [33] the spin of the atom 
is no longer constrained by the magnetic field and, due to 
the additional atomic spin degrees of freedom, the BEC 
exhibits a richer spinor order parameter structure. The 
spin of the optically trapped BECs can generally have sig- 
nificant effects on the dynamical properties of the BECs 
[33I , [sj, give rise to spin textures 135]^ an d simport of 
highly nontrivial defect structures [^, [s^, IH, [39[ . Ex- 
periments have also explored dynamics in spinor conden- 
sates in harmonic traps [10, [4l|, [42|, HI] and optical lat- 
tices [44] , and the application of spinor gases to spatially 
resolved magnetometry [45]. 

In this paper we investigate both dynamical and en- 
ergetic instabilities in the transport of multi-component 
BECs in optical lattices. We first consider magnetically 
trapped two-component BECs where the two conden- 
sates occupy different hyperfine states of the same atom 
or are formed by mixtures of two different atoms. Trans- 
port properties of two-component BECs in an optical lat- 
tice were studied in Ref. ^], and numerical results for 
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dynamical instabilities were presented for a special case. 
In particular, dynamical instabilities were shown to arise 
from the critical velocity as well as from the phase sep- 
aration of the two species. In contrast to that work, we 
obtain analytic results for the condensate dynamical in- 
stability points and analyze in detail the complete phase 
space of stability criteria for both dynamical and ener- 
getic instabilities across a broad range of parameters. We 
vary the intra- and inter-species interaction strengths, 
the site hopping term for each spin component indepen- 
dently, as well as the velocities of the two BECs, allowing 
application of our results to a large variety of experimen- 
tal systems. Among the novel results presented here are 
the possibility of a second BEC component stabilizing 
the superfluid flow of an otherwise unstable first BEC 
component (that exceeds the critical velocity of a single- 
component BEC). In addition, we find the free space 
phase separation criteria, that the square of the inter- 
species interaction coefficient exceeds the product of the 
intra-species interaction coefficients {U12 > ^11^22), can 
be reversed in an optical lattice. This can happen if 
one of the BECs has a velocity larger and the other one 
smaller than the single-species critical velocity (the ef- 
fective masses of the two components exhibit different 
signs). 

We also analyze transport properties of optically 
trapped spin-1 BECs in optical lattices, which have not 
been experimentally investigated to date. Here we ob- 
tain analytic expressions for both the dynamical and en- 
ergetic instability regions of the Bloch wave solutions. In 
contrast to the two-component case, spin changing colli- 
sions allow the atom population of different spin com- 
ponents to adjust to lower the energy of the system, 
according to whether the scattering lengths correspond 
to polar or ferromagnetic values. Our results illuminate 
the different stability properties of the polar versus fer- 
romagnetic solutions, which, in the absence of the Zee- 
man shifts, are most apparent for large spin-dependent 
scattering lengths or when the spin-dependent and spin- 
independent scattering lengths exhibit different signs. 
The presence of the Zeeman level shifts provides a richer 
variety of steady-state Bloch wave solutions, including 
novel solutions that do not exist for the case of small 
level shifts. We find that the quadratic Zeeman shift, 
due to its role in the energy conservation of spin chang- 
ing collisions, plays an important role in the stability of 
various condensate solutions. However, we also see that 
linear Zeeman shifts play an important role in stabiliz- 
ing many of the solutions. The dynamical instabilities 
of the spinor BECs can be important, e.g., also in the 
formation of solitons that have been studied in the ho- 
mogeneous case in Ref. (47| . 

In Sec. [ni we introduce the discrete nonlinear 
Schrodinger equation and the Bogoliubov-de Gennes ap- 
proach to study the stability properties of condensate so- 
lutions. This is done in the context of the two-component 
case but the same method is used later for the spinor case. 
We then derive expressions for the normal mode energies 



and discuss the dynamical and energetic stability of the 
two-component case. In Sec. Illll we apply this method to 
the spinor case. We first discuss the polar case, then the 
ferromagnetic case, then finally the effect of the Zeeman 
shifts. Some experimental considerations for observation 
of the effects studied are discussed in Sec. llVi We sum- 
marize our results in Sec. [Vl The Bogoliubov-de Gennes 
matrices are presented explicitly in Appendix El and the 
detailed analysis of the stability of the two-component 
case when the phase separation condition is reversed is 
given in Appendix [Bl 



II. A TWO-COMPONENT CONDENSATE IN 
AN OPTICAL LATTICE 

A. Two species system description 

Two-component BECs can be prepared in magnetic 
traps by simultaneously confining different atomic species 
in the same trap. The atoms may occupy two differ- 
ent hyperfine states of the same atomic species or form 
a mixture of two condensates of two different atomic 
species. For instance, two BEC components in perfectly 
overlapping isotropic magnetic trapping potentials were 
experimentally realized in hyperfine spin states of ^^Rb, 
I T) = |F = 2,m/ = 1) and | i) = \F = l,m/ = -1). 
In this system the inter- (a^i) and intraspecies (a^^ 
and ail) interaction strengths are nearly equal, with 
ail '• ' ^TT •* 1-024 : 1 : 0.973 [48]. Since the scattering 
lengths satisfy a^^ ^ anan, the two species experience 
dynamical phase separation and can strongly repel each 
other ^] . A more strongly repelling two-component sys- 
tem of different species was created using a ^^K-^^Rb 
mixture [26]. The interatomic interactions of two mag- 
netically trapped BEC components do not mix the atom 
population and the atom numbers of both species are 
separately conserved. 

The dynamics of the BECs follow from the coupled 
Gross-Pitaevskii equation (GPE) 

= (-^v^ + W + E ■ (1) 

k 

Here we have defined the interaction coefficients f^u = 
4:7rh^aii/mi and tZij = 2TT?i^aij/ii [i 7^ j), where the 
wavefunctions are normalized to Nj^ N = Ni -\- N2 is 
the total atom number, rrij is the atomic mass of BEC 
component and /i = mim2/(mi +7712) is the reduced 
mass. The intraspecies and the inter species scatter- 
ing lengths are denoted by an and a^j {i 7^ j), respec- 
tively. The external potential is generally a superposition 
of a harmonic trapping potential y^^^(r) = m{ujj^x'^ + 
^jyV'^'^^jz^'^)/'^ periodic optical lattice potential 

Vi^'\r) = V^'^ sm'{7Tx/a+^j), Vj{r) = vjj'\r) +V^^\r) , 
where a denotes the lattice spacing. In the following we 
ignore the effect of the harmonic trapping potential along 
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FIG. 1: The two BECs may experience different optical lattice 
potentials that can be shifted with respect to each other. The 
interspecies interaction strength U12 is proportional to the 
overlap integral of the two lattice site mode functions and 
it can be adjusted by shifting the lattices. Moreover, the 
kinetic energy hopping amplitude of the two species Ji , J2 
may be independently modified by changing the barrier height 
between the neighboring sites. 



the lattice and consider the system as translationally in- 
variant. We also neglect density fluctuations orthogonal 
to the optical lattice and consider the dynamics as effec- 
tively ID. 

We write the GPE in the tight-binding approximation 
by expanding the BEG wavefunctions on the basis of the 
Wannier functions and only keep the lowest vibrational 
states in each lattice site 77, so that '0j(r) = c^^ (j)jr^(Y) 
[49] . We obtain discrete nonlinear Schrodinger equations 
(DNLSEs): 



iri — 
dt 



(2) 



With similar assumptions, we have the hopping ampli- 
tude, Jj > 0, for the atoms between adjacent lattice sites: 



2m 



The nonlinearities 



are given by Ujk — 
f^jrji \^kr]\~- The two BEG species may 
generally experience different lattice potentials ^ and 
so the interspecies coupling coefficient U12 may be varied 
by displacing the two lattice potentials with respect to 
each other (to modify the overlap integral between the 
wavefunctions), even when the values of the scattering 
lengths remain constant; Fig. [H 



B. Dynamical stability of two species 

1. Collective two-component excitations 

We study the stability of plane wave solutions to 
Eq. ([2]) by investigating the effect of small perturbations 
around the carrier wave. Our treatment is analogous to 
the approach in Ref. [l^ to analyze single-component 
BEGs. For the constant atom density along the lat- 



tice, the Bloch waves Cr 



exp [iikjarj — i^jt)] that 



satisfy Eq. ([2]) exhibit the frequency = ^^riiUj 
2 Jj cos ( j a) , where rij 



.0-)|2 



denotes the constant 



atom population in spin state j at each site. The per- 
turbed carrier wave can be written as a Bogoliubov ex- 



pansion: 



(4) 



Here the kj represent the potentially non-zero velocities 
of the condensate in each component j. The stability of 
such states in a single component BEG has been studied 
experimentally in Refs. [?, 8] by dynamically moving the 
lattice potential. Substituting this expansion into Eq. ([2]) 
and linearizing the equations in lii, ^i, 2x2, and V2 yields 
a system of four equations 



dt 



\V2J 



a = 



a, 




(5) 



where (Jz denotes the 2x2 Pauli spin matrix. The el- 
ements of the 4x4 matrix A^(^) follow from the lin- 
earization procedure as in the single-component BEG 
case. We write out this matrix explicitly in Appendix 
[XI The eigenvalue problem for the matrix a^Aicj) can 
be solved analytically. In order to preserve the symme- 
try properties the Bogoliubov equations and to obtain 
simple analytic expressions for the normal mode energies 
we require that the atom currents of the two BEGs are 
equal, i.e., J\ sin(/cia) = J2 sin(/c2a). 

The eigenvalues represent the normal mode energies 
and read 



2Ji sin(/cia) sin((7a) ± W ^(^i^^ + ± \\J {'^i.q ~ ^2,qY + 16ei,q cos(/cia)e2,q cos(/c2tt)^i^2^i2 • (6) 



The four eigenvalues correspond to all permutations of ergies (without the Doppler shift term) 
the ± signs. Only two of the eigenvalues are independent. 
The first term in Eq. (|6]) represents the Doppler shift 
of the excitation energies due to the superfluid current. 
Here uoj^q denotes the single-condensate normal mode en- 



= ^j^q cos{kja)[ej^q cos{kja) + 2njUjj] , (7) 
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and 

e,,, = 4J,sin2(f) , (8) 

is the spectrum of an ideal, non- moving BEC. 

For the case of positive definite M{q) all the eigen- 
values ujq of aM{q) are real. In that case the physical 
solutions of the corresponding eigenvectors y exhibit pos- 
itive normalization yVy = 1 (the sign in the front of 
the first square root) and unphysical eigenvectors nega- 
tive normalization yVy = —1 (the ' — ' sign in the front 
of the first square root). The eigenvalues ujq with a non- 
vanishing imaginary part are associated with eigenvec- 
tors satisfying yVy = 0. The BEC system becomes 
dynamically unstable when the normal mode frequencies 
in Eq. (|6]) exhibit nonvanishing imaginary parts, indicat- 
ing perturbations that grow exponentially in time. Such 
modulational instabilities occur in a closed system due to 
the nonlinear dynamics and do not require energy dissi- 
pation. The rate at which the instability sets in depends 
on the magnitude of the imaginary part of the eigenfre- 
quency. 

I 



The long wavelength excitations are unstable when one 
of the solutions for the speed of sound has an imaginary 
part. 

If we do not assume that the two BEC currents are 
equal Ji sin(A:ia) 7^ J2 sin(/c2a), the Doppler shifts for the 
two BECs are different and the usual symmetry proper- 
ties between the positive and the negative energy Bo- 
goliubov eigenfunctions are lost. We still find analytic 
solutions for the eigenenergies uOq, but these no longer 
have simple compact expressions as in Eq. ([6]). The ba- 
sic formalism may be used to generate stability diagrams 
in these cases numerically, for instance, even if the two 
BECs have the velocities in the opposite directions. In 
the following we concentrate on analyzing the general fea- 
tures of the two-component system that may already be 
obtained from Eq. (|6]). 



2. Stability with equal signs for cos(/cia) and cos{k2a) 

We first analyze the dynamical stability of the two- 
component system, given in Eq. (|6]) for Ji sin(/cia) = 
J2sin(/c2a), for the case that cos(/cia) and cos(A:2a) ex- 
hibit equal sign. For that case, the expression inside the 
square root in Eq. (|6]) is always negative (for any values 
of and the dynamics unstable, if (jf^q + cj^^^ < 0. 

Physically, this corresponds to the situation where the 



For small momenta, qa <^ 1, cj^q ^ Jjq^o? = 
?i^g^/2m*, where we introduced the effective mass of a 
noninteracting BEC as m* = / {2JjO?). Similarly, we 

obtain 2Ji sin(/cia) s\ii{qa) :^ 2Jikiqa? = h^kiq/ml rein- 
forcing the interpretation of the first term in Eq. (|6]) as 
the Doppler shift contribution. 

If we set U12 = in Eq. (|6]), we obtain inde- 
pendently the decoupled normal mode energies of the 
two BECs ujq = 2Ji sin(/cia) sin(g'a) -j- \uji^q\ and ujq = 
2Ji sin(/cia) sin(ga) -h |^2,g|, analogously to the single- 
condensate normal modes obtained in Ref. [l^]- 

The intraspecies interaction U12 mixes the normal 
modes of the two BECs. In the experimentally interest- 
ing regime njl/jj ^ Jj, if 17^2 — ^11^22, one of the fre- 
quencies approaches zero indicating an instability similar 
to the uniform two-component BEC system. Specifically, 
for ki = k2 = 0, we obtain in that case uo'^^j^ :^ +^2,g 
and ul_ < ul^q.ulq. 

By expanding uOq for small q in Eq. (j6|) with ki = k2 = 
0, we obtain 00 q :^ hsq where s is the speed of sound 



(9) 
I 

two-component dynamical instability is driven by the in- 
stabilities of the individual single-component BEC exci- 
tations ([7j). We have cc;^^ + cc^^^^ < 0, for some values of 
q, if 

• 2 fQ^\ Dii + D22 .^^x 

V 2 y 2[J2 cos2 {kia) + J| cos2 {k2a)] ' ^ ' 

where 

Dij = JiCos{kia)njUjj . (11) 

The inequality is most easily satisfied for excitations 
in the long- wavelength limit <7 ^ 0, and is satisfied when 
the right hand side is positive, i.e., for 

^11+^22 < 0. (12) 

According to Eq. ([Ej), the modes can become unstable 
if kia,k2a > 7r/2 even when Uii,U22 > 0, which is the 
usual high velocity instability seen in the single compo- 
nent case [HI. The dynamics can also be unstable when 
kia,k2a < 7r/2, for negative values of Uu and U22- 

In addition to the instability occurring for ^i^q-\-^2,q < 
0, the two-component system in Eq. (|6]) is always dynam- 
ically unstable (for the case that cos(/cia) and cos(/c2a) 
exhibit equal sign) if 

U!2>UuU22, (13) 



s± = -\l JiUiUii + J2n2U22 ± V {JiniUii - J2n2U22Y + 4Ji J2nin2C/i2 • 
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as the expression in the outer square root in Eq. (|6]) be- 
comes negative. In particular, the unstable q modes are 



those that satisfy (when uo\^q + uo^^q > 0): 



J 



sm 



(f) 



< 



-{Di2 + D21) + \/ {D12 - D21Y + 4Ji cos(/cia) J2 cos{k2a)nin2Ui2 
4Ji cos(/cia) J2 cos(/c2<^) 

I 



(14) 



The two-component system, consequently, is in this case 
dynamically stable if Dn + D22 > and U12 ^ ^11^22- 

The unstable dynamics for U12 > ^11^22 correspond to 
the analogous instability which occurs in the free space 
case due to phase separation. One should emphasize, 
however, that the value of /7i2 is not only determined by 
the inter-species scattering length, but also by the spatial 
overlap integral of the lattice site wavefunctions for the 
two species; see the definition of Uij below Eq. Q. By 
means of shifting the relative position of the two BEG 
lattice potentials, one may easily reduce the value of /7i2. 

Phase space diagrams of the dynamical instability 
strengths as a function of the excitation wavelength q 
and condensate wave number k are shown in Fig. [2l In 
Fig.[2fa) we see the system is stable for the case ncUi2 = 
8J1, riiUii = 7121122 = 10 Ji (defining ric = y/nin2) when 
ka < O.Stt, then becomes unstable for ka > O.Stt, in ac- 
cordance with Eq. (p!2|) . The strength of the instabilities 
(the largest imaginary part of the eigenvalues) are linear 
in q [for the slope, see Eq. (|9])] until they saturate approx- 
imately at a value Im(cjg) c:^ \/ —^^ncU 12 J i cos (/cia), 
which is ~ llJi at /ca = TT in Fig [21(a). The strongest in- 
stability at /ca = TT represents period doubling that drives 
the system away from the Bloch state. The behavior 
is qualitatively similar for all nc/7i2 < tliUh = 77,2/722 
and also for smaller values of Ji and J2. However, for 
ncUi2 > niUii [Fig. EJb)], we see the predicted phase 
separation instability arises in the ka < O.Stt region. It 
is interesting to note this instability is markedly weaker 
than the high velocity instability {ka > O.Stt), as its 
maximum value (occurring for A: = 0, = tt) scales 
as lm{u;q) \/S{ncUi2 — niUii)Ji (in the limit that 
^c^i2 — ^1^11 ^ -^i)- In the case plotted in Fig. [2](b), the 
value saturates at ~ 3Ji, whereas the instability strength 
in the ka > O.Stt region reaches ~ 13Ji. 

Figures[3](a,b) compare the phase separation instability 
versus ncUi2 for kia = k2a = and kia = k2a = OAtt. 
Generally speaking, the larges t imaginary value re aches 
a maximum value lm{cOq) \/S{ncUi2 — niUii)Ji. The 
dependence on the sign of /7ii,/722 is shown in Fig. [H 
We see in Fig.Hl^a) that, for ncUi2 < riiUu = 712^/22 and 
cos(/cia) > 0, the instability is restricted to the attractive 
cases Uii < 0. Fig. (H^b) demonstrates this instability 
switches to the repulsive case for cos(/cia) < 0. Again 
these instabilities reach strengths lm{ujq) c:^ ^IGniUnJi. 
However, for the cos(A:ia) < case there occurs also a 
much weaker instability for attractive interactions with 
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FIG. 2: (a) Largest imaginary part of the eigenvalues Eq. JB]) 
as function k = ki = k2 and q for the case niUn = n2t/22, 
ncUi2 = 8J1 and J2 = Ji. We define ric = ^n\n2. Gray 
scale goes from (black) to 5Ji (white), (b) Same plot with 
ncUi2 = 13Ji (phase separation regime). 



strength Im(cjg) :^ ^y8{ncUl2 — niUii)Ji. 



3. Stability with 



signs for cos{kia) and cos(k2a) 



The case that cos kia and cos/c2a have different signs 
represents a configuration where the velocities of the two 
BEGs are located on the opposite sides of the defiection 
point in the ideal, single-particle BEG excitation spec- 
trum ([8|) (the effective masses of the two components 
exhibit different signs) and is presented in detail in Ap- 
pendix [B1 In that situation, we also always find a dy- 
namical instability when ^2,q < resulting in 
the relation similar to Eq. p2|) . In this case, however, 
the high velocity instability condition is highly nontriv- 
ial, depending on the values of the hopping amplitudes, 
the interaction strengths, atom numbers, and the veloci- 
ties: One of the BEGs that reaches the single-component 
critical velocity kia > 7r/2 may, or may not, destabi- 
lize the two-component BEG system, depending on the 
parameter values. Moreover, for cc;^^ + cj^^g > 0, the dy- 
namical stability condition due to the phase separation 
of the non- moving system, U12 < ^11^22, is reversed, so 
that the entire dynamically stable region occurs, when 
U12 > U11U22' In particular, we find in that case the 
system to be dynamically stable if U12 satisfies, depend- 
ing on the value of t/n, either U11U22 < ^12 < oi" 
^11^22 < ^2 < U12 < ^1, where ^1 and ^2 are defined 
in Eqs. ()B3p and ()B8p . Interestingly, this also repre- 
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FIG. 3: (a)-(b) Imaginary part of the eigenvalues versus 
ncUi2 (again ric = a/tiitt^), holding niUn = n2U22 = 10 Ji 
for the indicated condensate velocities k. Here J2 = Ji (c) 
A case with opposite signs of cos(/cia) and cos(/c2a), showing 
a region of dynamical stability for all q. Here niUn = 10 Ji, 
n2U22 = 12Ji and J2 = 2Ji. A similar narrow stable region 
exists around nct/12 — — 12Ji. (d) Another case with op- 
posite signs of cos(A:ia) and cos(A:2a). Here niUn = 27 Ji, 
^2^/22 = 30 Ji and J2 = Ji. Again another stable region in 
this case is located close to ncUi2 — — 28.5Ji. In these plots 
the gray scale goes from to 5Ji. 
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FIG. 4: Imaginary part of the eigenvalues versus q and varying 
niUii — n2U22, keeping rict/i 2 — niUn — Ji and with J2 — Ji. 
In this figure, the gray scale runs from to 2.5 Ji. 



sents a situation where the other condensate component 
can stabilize the superfluid flow of an otherwise unsta- 
ble condensate (exceeding the single-component critical 
velocity). Moreover, the two-component system may be 
dynamically stable even for Un < and U22 > 0, since 
in that case U11U22 < ^12- 

Figures[3Kc-d) show cases with cos(/cia) and cos(/c2tt) of 
different sign [but satisfying the condition Ji sin(/cia) = 
J2 sin(/c2a)]. For the parameters of Fig. ^c) there is a 



range of ncUi2 for which the system is dynamically sta- 
ble. In this case niUn satisfies Eq. ()B6p and the dy- 
namically stable region, according to Eq. (jBTj), is ^2 < 
U12 < £,1' The relevant quantities are nc\fQi — 11.5Ji 
and nc^/Qv = 12.5Ji for this example. Note that the 
dynamical instabilities tend to be much weaker on the 
large U12 side of this stability range. In Fig. [3jc) we 
have different hopping amplitudes Ji 7^ J2. Fig. [31(d) 
shows a case with a smaller range of stable ncUi2 but 
with Ji = J2. In this case niUn satisfies Eq. (|B4p and 
the dynamically stable region, according to Eq. ()B5p . is 
for U11U22 < U^2 < 6- Here ncVUiiU22 = 28.46 Ji and 
ncVTi = 28.5OJ1. 



C. Energetic stability 

The energetic stability of the superfluid flow of the ho- 
mogeneous two-component mixture depends on the prop- 
erties of the energy functional. The second-order varia- 
tions of the energy for small perturbations in the carrier 
wave are determined by the matrix Ai{q) and the system 
is energetically stable if M{q) is positive definite. If any 
of the eigenvalues of M{q) are negative, the system may 
relax to a state with lower energy by means of dissipa- 
tive coupling to the environment. The rate at which such 
relaxation happens depends on the strength of the cou- 
pling, e.g., on the number of thermal atoms interacting 
with the condensate. 

The eigenvalues of A4 {q) can also be evaluated analyt- 
ically but the full solutions are rather lengthy. In Fig. [5] 
we show the energetically unstable regions of the two- 
component dynamics. Note that positive definite M{q) 
implies real eigenvalues cOq of (jA4{q) in Eq. (|6|), so the 
dynamically unstable region always forms a subset of the 
energetically unstable region. Figures [5l(a-b) show cases 
with ncUi2 < niUii and ncUi2 > niUu^ respectively. In 
the latter case there are regions of instability for some 
q at all condensate wavenumbers k. In the former case, 
there are is a band of k with all q modes energetically sta- 
ble, with a width proportional to y^2Ji{niUii — nct/12), 
which determines the speed of sound for the spin wave; 
see Eq. (j9j) with riiUu = n2/722, Ji = J2- This is anal- 
ogous to the single-component case [12] where there is a 
band of energetically stable k with a width proportional 
to the speed of sound ex ^/2JiniUii. In Figs.[5](c),(d) we 
show the dependence on U12. It is seen that at /c = 
instability only occurs for U12 > Uu while for finite k 
the condition for stability becomes more stringent. For 
ka > O.Stt the entire region is unstable. We also cal- 
culated that the specific parameter regimes correspond- 
ing to dynamical stability with two different condensate 
wavenumbers /ci,/c2 (Figs. [3jc,d)) are energetically un- 
stable at all q. 
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FIG. 5: (a)-(b) Regions of energetic stability (white) and 
instability (black) versus ki and q for the cases indicated. We 
choose niUii = n2U22 = 10 Ji, J2 = Ji and k2 = ki. (c)-(d) 
Energetic stability regions versus ric t/12 for a stationary and 
moving condensate case. 



III. A SPIN-1 CONDENSATE IN AN OPTICAL 
LATTICE 

A. Spinor Gross-Pitaevskii equations 

We now consider a BEG of spin-1 atoms. In the ab- 
sence of a magnetic trapping potential, the macroscopic 
BEG wave function is determined by a spinor wave func- 
tion ^ with three complex components |50]- The Hamil- 
tonian density of the classical GP mean-field theory for 
this system reads: 



«=^|V*|» 



■5i{B-F)/j + 52{(B 



C2tr 
2 



(15) 



In Eq. (p!5|) , F is the vector formed by the three com- 
ponents of the 3x3 Pauli spin-1 matrices [50|, (F) = 
. F . denotes the average spin, and p(r) = |^(r)p 
the total atom density. The weak external magnetic field 
is denoted by B and is assumed to point along the z axis. 
The magnetic field produces the linear and quadratic Zee- 
man level shifts whose effect is described by the last two 
terms in Eq. ([15]). As in the two-component case, the ex- 
ternal potential V is the sum of the harmonic part Vh (in 
this case due to the optical dipole trap) and the optical 
lattice potential Vl: V{r) = l/;/(r)+14(r). In the follow- 
ing we ignore the harmonic potential and, for simplicity, 
assume that the lattice potential is the same for all the 
spinor components. Thus the Wannier basis functions 



(pr] of our discrete basis no longer depend on the internal 
state j. Here cq and C2 are the spin-independent and spin- 
dependent two-body interaction coefficients. In terms of 
the s-wave scattering lengths ao and a2, for the channels 
with total angular momentum zero and two, they are: 
Co = 47rfi^(2a2 + ao)/3m, and C2 = 47rfi^(a2 — ao)/3m. 
For ^^Na, (a2 - ao)/3 0^ 2aB and (2a2 + ao)/3 SOa^, 
where as = 0.0529 nm is the Bohr radius [50], indicating 
C2/C0 — 0.04. In contrast to the two-component case, in 
which the atom number in the two-components do not 
mix, in the spinor case a homogenous condensate wave- 
function can adjust itself by varying the relative atom 
populations. From Eq. ([15]) in the absence of the ex- 
ternal magnetic field we immediately observe that, since 
C2 > for ^^Na, corresponding to the polar phase, the 
energy is minimized by setting (F) = throughout the 
BEG for the case of a uniform order parameter field. Al- 
ternatively, for ^^Rb we have C2/C0 ^ -0.0036 [si^. The 
parameter values for ^^Rb correspond to the ferromag- 
netic phase, since C2 < 0, and the energy in Eq. ([T5|) in 
the absence of the external magnetic field is minimized 
when 1(F) I = 1 throughout the BEG for the case of a 
uniform spin distribution. 

Again using the lowest band of the Wannier state basis, 
the DNLSEs are written: 



ih 



dt 



(+) 



..{+) 



Uo 



E 

a=+,0,- 



.(-)|2 



.(0)|2>) (+) 



77 



j(4;U4-i) + '^-4"^ 



a=+,0,- 



'r]-lJ 



T] I 

* (0)2 
77 ^77 



.(+)|2 



.(0)|2>, (-) 

-"77 I 7^77 



■U2ci+^-c: 



ih 



(0) 



dt 



^+1 



Q! = + ,0,- 
+ t^2(|cWr 



..(-)\2\J0) 



V 



(16) 



Here J is defined as before ([3]), but with no de- 
pendence on the internal state j [HH and /7o,2 — 
Co, 2 / \4^'n\^' The primary qualitative difference with 
the two-component case, seen in the last term of each of 
these equations, is the allowance of spin exchange colli- 
sions. The additional energy shift terms ^+,^- account 
for Zeeman shifts (with respect to the level m^? = 0) due 
to an external magnetic field B. For simplicity, we ignore 
any effects of magnetic field gradients. 
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We study the stability of moving Bloch wave solutions 
to the DNLSEs (p!6|) . In order to find the low energy 
stationary solutions, we substitute 



nexp [z(/ca77 — /it)] . (17) 



expressions for the normal mode energies 





Here /i is the chemical potential and n = X1q;=+ o - k^*^^ P 
is the total condensate density that is assumed to be 
constant along the lattice. The spinor wave function, 
= (C+ , Co 5 C- ) 7 satisfies the normalization condition 

■ C = 1. We concentrate on solutions for which ( is 
constant along the lattice. 

We substitute the same Bogoliubov expansion as in 
Eq. (HI) for the linearized fiuctuations around the carrier 
wave solution (pT|) in the DNLSEs ([16]). This yields a 
6x6 matrix aM, analogous to Eq. ([5|) (in this case with 
a having three Pauli matrices in the diagonal), gov- 
erning the dynamical stability of the system. The matrix 
Ai is given explicitly by Eq. ()A2p in Appendix [Al 

As in the two-component BEC case, negative eigen- 
values of the matrix M. indicate the regions of energetic 
instability, while the eigenvalues of aM yield the normal 
mode frequencies. The imaginary parts of these normal 
mode frequencies represent the strength of dynamical in- 
stabilities. 



B. Stability in the polar case 

1. Dynamical stability 

In the polar case (that is energetically favored for U2 > 
0), we consider uniform spin profiles with the average 
spin value zero, and assume no Zeeman shifts for the 
time being (5+ = S- =0. All the degenerate, physically 
distinguishable, ground state solutions for U2 > may 
then be determined by means of the macroscopic BEC 
phase cp and a real unit vector d defining the quantization 
axis of the spin. The spinor wavefunction reads: 




(18) 



As in the similar polar phase of superfluid ^He-A [53|, 
the states (d, cp) and (— d, cp -\- tt) are identical. This can 
be conveniently taken into account by considering the d 
field to define unoriented axes rather than vectors. 

The solution (pTj) with Eq. (p!8|) has a chemical poten- 
tial value fi = —2Jcos{ka)-\-Uon for any choice of (d, (p). 
Since also the excitations are the same for any values of 
(d, (/?), we may choose the simplest form of the matrix 
M in Eq. ()A2p . that is obtained by choosing d to point 
along the z axis and cp = 0. 

By calculating the eigenvalues of aM. with this partic- 
ular choice of the BEC wavefunction we obtain analytic 



ui±{q) = Cq^k ^ y cos{ka) [cq cos{ka) + 2nUo] 



^2±{q) = Cq^k ^ y cos(/ca) [cq cos(/ca) + 2nU2] (19) 

where UJ2± are each doubly degenerate. The physical 
solutions correspond to the sign in the front of the 
square root. Here again Cq denotes the spectrum of an 
ideal, non-moving BEC 

e,=4Jsm'{^) , (20) 

and the Doppler shift term in the energy is given by 



Cq^k — 2 Jsin((7a)sin(/ca) 



(21) 



It is clear from Eq. (p!9|) that, for U2 < Uq and Uq > 
0, uji-^{q) drives the instability. The uji-^{q) modes are 
unstable when the expression inside the square root is 
negative, which happens for q values that satisfy: 



sm 



(?) 



< 



nUo 



2 J cos {ka) 



(22) 



At least some modes are unstable whenever ka > 7r/2 
and all the q modes are unstable when —nUo/{2J) < 
cos {ka) < 0. Figure [6ja) plots the largest imaginary 
part of the mode frequencies Eq. ([19]) versus q and k for 
the case nUo/J = 100 and nU2/J = 4 (corresponding 
to ^^Na). As in the two-component case, one sees the 
system is stable for ka < 7r/2, while for ka > 7r/2 one 
has an instability with a growth rate linear in q in the 
long wavelength limit before saturating at a maximum 
value ^ ^y—8nUoJ cos {ka) ^ 27 J. Figure [61(b) plots a 
case with a much smaller nonlinearity {nUo/J = 1 and 
nU2/J = 0.04). While ka < 7r/2 is still the condition for 
stability of all the modes, one sees that for higher values 
of k there exists only a band of low unstable modes in 
the lower q region. This is due to the fact that the RHS 
of ([22]) becomes less than unity and thus can be exceeded 
by the LHS for large q. 

Just as we saw in the two-component case (see Fig. [4]), 
these conditions are somewhat reversed for attractive in- 
teractions Uq < 0, as then some excitations of u;i-^{q) are 
unstable whenever ka < 7r/2 and all the q modes are un- 
stable when < cos {ka) < —nUo/{2J). The dependence 
on the interaction coefficient Uq is shown in Fig. [3 In 
these plots, we kept U2n/J = 2 constant. For ka = 0.27r 
(Fig. [TJa)) unstable modes occur for negative /7o, while 
for ka = O.Ttt (Fig. [71(b)), this instability occurs for repul- 
sive interactions Uq > 0. Also, as in the two-component 
case, there is an additional, weaker instability in the at- 
tractive case Uo < with ka > 7r/2. This instability 
is driven by uj2-\-{q) and has the same condition for in- 
stability ([22]) with Uq replaced by U2. The magnitude 
generally reaches ~ \^SnU2J for nU2 ^ J- In the case 
plotted in Fig.[71(b), it has a maximum magnitude ^ 2 J. 
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FIG. 7: Imaginary part of the eigenvalues versus n[/o, for 
nU2 = 2 J for the condensate velocities indicated. Gray scale 
goes from to 5 J. 



FIG. 8: (a) Energetic stability (white) and instability (black) 
regions for tiUq = lOOJ, nU2 = 4J. (b) Dependence of 
the energetic stability of interaction coefficient n[/o, keeping 
nU2 — 2J. The behavior remains qualitatively the same for 
non-zero ka < 7r/2, but with the instability region extending 
into the ?7o > region. 



modes for cos(/ca)/ sin^(/ca) < 2J/nU2' For small veloc- 
ities, ka <C 1, this c ondition is approximately equal to 
Tik/m'' > Y^n/72/m*, where m* = Ti^ /(2Jo?) is the effec- 
tive mass of a noninteracting BEG. This energetic insta- 
bility threshold demonstrates the Landau criterion that 
the velocity becomes larger than the speed of sound (of 
spin waves) a^/2JnU2/Ti = ^/nU2/m'^. For U2 > Uq the 
instability driven by ei_ and the condition is the same, 
but replacing U2 Uq. When Uq becomes negative, 
there is an additional instability, as shown in Fig. [8](b). 



The equivalence of the mode energy dependence on 
U2 and Uo in Eq. (p!9|) implies that whenever Uo and U2 
are of opposite sign and much larger than J, at least 
one of the eigenvalues will be imaginary at any k. In 
the previous paragraph we discussed how this resulted 
in an instability for attractive condensates Uq < with 
polar spin-dependent scattering lengths U2 > 0. Another 
implication of this is that polar-like condensate solutions 
(p!8|) with Uo > are unstable for ferromagnetic spin- 
dependent scattering lengths U2 < 0. Only when both 
are positive or both are negative is there a region of k 
with dynamical stability. 

2. Energetic stability 

Turning now to the energetic instabilities, we also ob- 
tain analytic results for the eigevalues of M. and look for 
negative eigenvalues. We find 



^2±{q) = eq cos{ka) + nU2 ± ^C^^ + n2/7|. (23) 

For U2 < Uo 62- drives the instability. Figure [8] plots 
the regions of energetic instability. There exist unstable 



C. Stability in the ferromagnetic case 

1. Dynamical stability 

In the ferromagnetic case (that is energetically favored 
for U2 < 0) we consider uniform spin profiles for which 
the magnitude of the spin is maximized, |(F)| = 1, which 
minimizes the mean field energy (fT5]) . We again assume 
no magnetic Zeeman shifts J+ = (5_ = 0. As in analogous 
states for superfluid liquid helium-3 [H^ , the rotations of 
the spinor axes can be used to couple physically distin- 
guishable ground states. Here all the degenerate states 
are related by spatial rotations of the atomic spin axes 
and we may parametrize the spin wavefunction as 

C+\ /e— cos2(/?/2)\ 

Co =e^M sin(/?)/72 , (24) 
C-J \e'''sm^{/3/2) J 

where a, /3, (j) are the Euler angles. The solution (pT|) with 
Eq. has a chemical potential /i = — 2 Jcos(/ca) + (/7o + 
U2)n for any chosen ground state in Eq. ([24|) . The sim- 
plest form of the Bogoliubov-de Gennes matrix aM. may 
be obtained by choosing /3 = a = (j) = and substituting 
Eq. ^ into Eq. (|A2l) in Appendix [Al 
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The mode energies are found to be: 



^i±(g) = C'q,/c ± eg cos{ka) 
^2±{q) = Cq^k ± eqCos{ka) ^ 2nU2 
^3±{q) = Cq^k 



±Jeq cos{ka) [cq cos{ka) + 2n{Uo + ^2)] (25) 
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(c) Energetic stability versus the interaction coefficient Uo, 
keeping U2 = —2 J. 



2. Energetic stability 



The dynamical instabilities are driven entirely by cos-^ 
and the only difference with the polar case Eq. ([19]) is the 
replacement Uq and U2 individually by the sum t/o + ^2- 
This dependence can be understood from the fact that 
with 1(F) I = 1 the total nonlinearity is oc /7o + U2 and so 
this quantity determines the attractive or repulsive char- 
acter of the condensate. Thus for Uq > and \ U2\ <C |/7o|, 
the instability diagrams qualitatively similar to Fig. [6l 

Differences between the ferromagnetic and polar cases 
become clear when one examines the instability depen- 
dence on Uq. This is shown for the ferromagnetic case 
in Fig. [9l where we keep nU2 = — 2J constant, and 
should be contrasted with the polar case, Fig. [3 For 
ka < 7r/2 (the figure shows ka = 0.27r), an instability oc- 
curs for Uo < —U2 and increases with greater \Uo\ while 
for ka > 7t/2 (the figure shows ka = O.Ttt), the insta- 
bility occurs for Uq > — C/2. An important difference 
from the polar case is that, for ka > 7r/2, there is no 
instability for Uq < 0. In addition, the instability bor- 
der occurs at Uq = —U2 rather than at Uq = 0. Finally, 
from Eq. ([25]). we note that a ferromagnetic BEC solution 
([24|) with Uo > and a polar spin-dependent scattering 
length U2 > can be dynamically stable, in contrast to 
a polar solution with ferromagnetic scattering length, as 
discussed above. 



The energy eigenvalues in the ferromagnetic case are 

ei±(<?) = eqCos{ka) ±Cq^k 
^2±{q) = eqCos{ka) - 2nU2±Cq^k 
^s±{q) = Cq cos{ka) ^ n{Uo ^ U2) 

±^Cl,^n^{Uo^U2y. (26) 

Unlike the polar case, there is one energy eigenvalue 
ei_ corresponding to a pure (Doppler-shifted) kinetic 
energy. This gives rise to energetic instabilities for 
sin^(g'a/2)/ sin(g'a) < tan(/ca), as plotted in Fig. [TOlf a). 
Though no dynamic instability exists except for much 
higher /c, ferromagnetic spinor BECs are subject to this 
energetic instability in the presence of thermal excita- 
tion for any non-zero k. We note that for ka > O.Stt all 
q modes are unstable. 

Finally, for attractive condensates {Uq -\-U2 < 0) there 
are additional regions of energetic instability from 63-. 
Fig. [TOlfb-c) shows the dependence versus Uq. For k = 
there are unstable modes for all attractive condensate 
scattering length cases. A moving condensate {ka = 0.27r 
is shown in Fig. [TOlf c)). increases the region of unstable q 
modes for attractive condensates. The band of energetic 
instability at low q for /7o > in Fig. [TOlfc) is simply the 
Doppler induced energetic instability discussed above. 
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D. Effects of Zeeman splitting 

When Zeeman splittings due an external magnetic field 
[S± in Eq. (p!6|) ] are non-zero, the symmetry of the polar 
and ferromagnetic solutions, Eqs. (p!8|) and ([24|) , breaks 
down and we find a new set of steady-state solutions. 
Here we examine these solutions and again calculate 
the dynamic stability of these various solutions, particu- 
larly noting how the stability varies with the linear and 
quadratic Zeeman shifts, which we denote, respectively, 
as ^ = ((5+ - S-)/2 and J = ((5+ + (5_)/2. The linear 
Zeeman shifts are S = 2gFlJiBB = 2((27r)1.4 MHz/G)5, 
where gp is the Lande factor, and is —1/2 for the ground- 
state F = 1 manifold of ^^Rb and ^^Na. For Zeeman 
shifts substantially smaller than the hyperfine splitting, 
which is our interest here, the quadratic shifts S are typ- 
ically smaller than the linear shifts and can be extracted 
from the Breit-Rabi formula (sj]. For alkali atoms the 
quadratic shift is positive, but it can generally be of ei- 
ther sign. The level shifts in a spin-1 BEG may also be 
engineered in other ways, e.g., by using off-resonant mi- 
crowave fields that generate electromagnetically-induced 
level splittings [HH , allowing essentially arbitrary experi- 
mentally prepared level shifts for 5 and 5. Note that the 
linear Zeeman shift does not affect the energy conserva- 
tion of a spin-changing collision \mF = 0,mF = 0) ^ 
Imp = +I,mi7 = —1), while the quadratic Zeeman shift 
does, and thus plays an important role in the stability 
properties. 

One of the steady-state solutions in the presence of the 
Zeeman splitting has the chemical potential /i = /io = 
— 2Jcos {ka) + nUo and reads: 



(27) 



This solution forms a subset of the polar solutions (p!8|) in 
the absence of the magnetic splitting, with the d pointing 
along the z direction. 

The stability is again analyzed by substituting the 
steady- state solution [Eq. ([27|)] into Eq. (|A2]) . By cal- 
culating the eigenvalues of aAi we obtain analytic ex- 
pressions for the normal mode energies, as in Eq. (p!9|) . 
Here uji±{q) remains unchanged in the presence of the 
Zeeman splitting and 

Pq,k = V [Ca COS 



{ka) + 5] [cq cos(te) + 2nU2 + ^] , (28) 



The linear splitting lifts the degeneracy between 002±{(l) 
and CJ3± {q) in Eq. (p!9|) and the quadratic splitting intro- 
duces an energy gap 5 in the single-particle phonon mode 
spectrum cos(/ca) in ij02±{q) and uo^±{q). 

The dynamical stability will then be governed by the 
sign of the square root argument of Pq^k and is seen to be 
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FIG. 11: Largest imaginary parts of the normal mode energies 
for the polar- like solution Eq. (|3Q)) . In all plots, Uo = 100 J 
and U2 = 0.04[/o (corresponding to ^^Na). Plots are on a gray 
scale of to 3 J. (a) Dynamic instability versus 6 for linear 
Zeeman shift S — 2J and k — 0. (b-c) Versus linear shift 6 
for quadratic Zeeman shifts 6 = 0.5 J and 8 J, with /c = 0. (d) 
Versus S with ka = O.SStt. 



unaffected by the linear Zeeman shift. The condition for 
the existence of unstable modes is —Cq cos{ka) — 2nU2 < 
S < —eqCos{ka) for U2 > 0, and —eqCos{ka) < S < 
—eqCos{ka) — 2nU2 for U2 < 0. If U2 > and ka < 
7r/2, Eq. ([28|) predicts dynamical stability for positive 
quadratic shift ^ > 0, and instability at some for ^ < 0. 
As in previous cases, we find much larger instabilities 
when ka > 7r/2 (for Uq > 0) from uji±. 

Similarly, we obtain the eigenvalues of Ai as in 
Eq. ([23|) . Here €i±{q) is unchanged and 



e2±{q) = Cq cos{ka) + nU2 -\- S ± 

^3±{q) = eq cos(ka) + n/72 + 5 ± \J {Cq^k - S)'^ ^ n^U^ . 

(29) 

We find another steady-state solution to the DNLSEs 
(fT6|) with the chemical potential /i = /io + ^ that reads 



e-'^\/l-S/nU2 




(30) 



■S/nU2 J 



This solution only exists for sufficiently small linear Zee- 
man shifts 1^1 < \nU2\. For ^ = 0, Eq. (|3Q|) coincides with 
the subset of solutions to Eq. ([18]). Although for small 
5 7^ the solution (|3Q|) is still close to the polar state 
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of Eq. ([18]), the spin expectation value is no longer zero, 
(F) = — (5/(n/72)z, with d restricted on the xy plane and 
the induced (F) pointing along the z axis [56|. At the 
boundary of the validity of Eq. ([T8|), 5 = ±n[/2, we obtain 



(F) = ^z, as in a ferromagnetic state, so that Eq. (|30|) in 
fact interpolates between the polar and the ferromagnetic 
solutions. Moreover, we again find analytic solutions for 
the normal mode energies: 



^i±(<7) 


= Cq^k ± \l 


' eq cos{ka) | 


Cq cos{ka) + 2n/72 — 


26] -] 






= Cq^k ± \J 


/ Cq cos{ka) 1 


Cq cos{ka) + n(/7o + 


U2)] 






= Cq^k ± \J 


/ Cq cos{ka) 1 


Cq cos{ka) + n([/o + 


U2)] 


~ Lq^k 




= Jel cos2 






]• 





(31) 



r 



For Uq^U2 > and /ca < 7r/2, whenever the solution 
(|3Q|) exists (i.e., when |n/72| > the dynamical in- 
stabilities are solely driven by u;i±{q). Moreover, under 
these conditions the mode uji±{q) is dynamically stable 
when Eq. ([3Q|) is energetically favorable to Eq. ([27|) (i.e., 
when S < 0). In Fig.[TT]we show some stability diagrams 
for the parameters of ^^Na. The instability dependence 
on the quadratic Zeeman shift for a particular linear Zee- 
man shift is shown in Fig. [TTlfa). This diagram may be 
understood by noting that the mode uji±{q) exhibits a 
nonvanishing imaginary part, if 

c_ < 5 < c+. 



c± = Cq cos{ka) + nU2 ± \J n'^U^ - S'^ . (32) 

This forms an instability stripe with a width 2(71^/71 — 
j2ji/2 Fig. [TTlf a) and for larger values of S the system 
stabilizes again. We observe the stripe position shift in 5 
by 4J from the qa = to qa = tt. We also calculated the 
energetic stability of Eq. ([30]) and found that in Fig.fTTTa) 
the region to the left of the stripe is energetically stable, 
while the entire region to the right of the stripe (i.e., large 
6) is energetically unstable. 

While not obvious in Fig. [TTlfa) there is a small region 
of stability for positive S. This is seen more easily by 
plotting the instability versus the linear Zeeman shift. 
In Figs. [TlTb-d) we plot this for the entire range of va- 
lidity of Eq. daO]) (1^1 < \nU2\) for various S. For small 
quadratic Zeeman shift, there exists a range of linear Zee- 
man shifts which stabilize the system, as it approaches 
the ferromagnetic state. For larger quadratic shifts this 
range shrinks, until eventually the system is unstable at 
all possible (5, as in Fig. [TTlfb). 

The effect of larger k is generally simply to stretch the 
region of instability to larger ranges of q for each case, 
as the effect of the Cq shift in Eq. (|32]) vanishes. For 
ka > O.Stt (with Uq > 0) the usual, and much larger, 
instability for large /c, seen in previous cases, dominates 
the stability diagram. Such a case is plotted in Fig. fTTT d). 




In the absence of the Zeeman splitting the polar state 
(p!8|) is always dynamically unstable for U2 < 0. For the 
state (|3Q|) . even close to the polar state, this is no longer 
the case. For Uq > and U2 < 0, Eq. (|32]) represents 
the entire unstable region, provided that 4^^|[/o//72| < 
n\Uo-U2)^ and \Uo\ > \U2l 

In the presence of the Zeeman splitting the ferromag- 
netic state is modified to 



(33) 



with the chemical potential /i = /io + + ^^+, or to an 
analogous ferromagnetic state with (("+, ^+) interchanged 
with ((^_,(5_). Although Eq. ([33]) is a steady-state solu- 
tion to the DNLSEs ([16]) for any values of (5, we also find 
that the solution ([3Q|) has the limit Eq. ([33]) at the bound- 
ary of the validity region 6 = —nU2- At the other limit of 
the validity of Eq. (|3Q|) (at 5 = nU2) we recover the other 
ferromagnetic state, defined by ((^_,(5_). Moreover, the 
solution (|33|) is energetically favorable to Eq. ([27|) when 
n/72/2 + (5+ < and to Eq. ([30]) when n/72/2 + ^ < 0. 

We find that the normal mode energies and the eigen- 
values of Ai corresponding to Eq. ([33]) are obtained 
from the non- Zeeman shifted mode frequencies, Eqs. ([25]) 
and ([26]) . by shifting the single-particle excitation ener- 
gies: eqCOs{ka) eqCOs{ka) — (5+ in uji±{q) and ei±{q); 
and eqCOs{ka) eqCOs{ka) — 2S in u;2±{q) and e2±{q). 
The energies u;s±{q) and es±{q) are unchanged. Because 
it is ujs± which drives the dynamical instability, the sta- 
bility diagram is unchanged from that of Eq. ([25]) . 

We find an additional ferromagnetic-like steady-state 
solution to the DNLSEs ([16]) with /i = /io + ^^2 + (^^ - 
S'^)/2S, that reads 



l2nU2S^P-P 



8nU2S^ 



Co = e^(7++7-)/2^i_|^^|2_|^_|2. 



(34) 
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FIG. 12: Largest imaginary parts of the normal mode energies 
for the solution Eq. (|34)) . In all plots, nUo — 100 J. Hatched 
areas represent regions where the condensate solution is not 
valid, (a)-(c) Dynamic instability versus linear Zeeman shift 
and q for U2 — — 0.0036[/o (corresponding to ^^Rb). Gray 
scale runs from (a) to 0.3J, (b) to 0.3J, (c) to 0.2 J. In 
(a) we emphasize that there is a range of ^, on each side of 
the instability region, stable for all q. (d) Dynamic instability 
versus spin-dependent scattering length [/2 for ^ = 2J and 
6 — J. Gray scale from to 0.2 J. 



This solution only exists if the expressions inside the 
square roots of Eq. ([34|) are positive, i.e., 



< 



2nU2^ + ^'^ 



SnU2S^ 



< 



1 



(35) 



This condition gives rise to a ranges of validity for U2 
in terms of two parameters = (^^ — P)/2S and 
/3+ = (^^ + ^^)/2^. In the most common case that the lin- 
ear Zeeman shift is larger in magnitude (|^| > \6\) these 
inequalities give a finite range: /3_ < n/72 < /^+ for J > 
and (3^ < nU2 < P- for J < 0. 

If 1^1 < the inequalities change direction, giving 
instead an intermediate range oi U2 where is not 
valid. In particular, the requirement for validity for ^ > 
is nU2 < P- or nU2 > while for ^ < it is nU2 < 
or nU2 > P-' 

At the lower limit of Eq. (|35]) {nU2 = P-) the solution 
coincides with the polar solution ([27|) with (F) = 0, 
while at upper limit {nU2 = /?+) it equals Eq. ([3Q|) . In 
general, the spin for the solution is non- vanishing 



l(F)P = 



5"^ - (5-* 



452n2j72 



(36) 



This condensate solution ([34j) is energetically favorable 
to (|33|) for positive quadratic shifts 5 > 0, the opposite 
relationship of solution (|3Q|) to ([27j) . 

We computed the dynamical and energetic stability 
of this solution. Figures [T2lf a)-(c) show the dynamical 
instability strengths as a function linear Zeeman shift 
for several quadratic shifts and Rb-87 scattering length 
{U2 = — 0.0036/7o)- The hatched areas and all \S\ larger 
than the range of these plots are regions where the con- 
densate solution (|34]) is not valid. For S > (Fig. [TSKa)) 
we see the valid regions where the solution exists are dy- 
namically stable. For ^ < (Figs. [T2lf b)-(c)) there are 
always some unstable modes q. However note the inter- 
esting behavior that for large \S\^ there exist only small 
bands of unstable q and weak instability (note the range 
of the plots in the caption). The system is energetically 
unstable for q in the regions overlapping and below the 
small dynamical instability bands. 

In Fig. [T2lf d) we plot the instability strength versus U2 
over its range of validity, which in this case is 1.5 J < 
nU2 < 2.5 J. Here we see a weak band of instability for 
low q in the range of U2 where the solution exists. 



IV. EXPERIMENTAL CONSIDERATIONS 

In the experimental realizations of optical lattice sys- 
tems, ultra-cold atoms have been trapped in a combined 
optical lattice and a harmonic trap. The transport prop- 
erties may then be studied by suddenly displacing the 
harmonic trap, e.g., by using a magnetic field gradient. 
This excites dipolar oscillations of atoms along the lat- 
tice direction with the maximum velocity proportional 
to the harmonic trap displacement [i, 0] . The other al- 
ternative is to use a moving-standing wave, so that the 
atoms are trapped close to the harmonic trap minimum 
and experience a moving optical lattice potential 0, @] • 
The advantage of the latter technique is that the velocity 
of the atoms with respect to the lattice is constant. In 
such transport experiments the dynamical instabilities 
may typically be observed on much shorter time scales 
than the energetic ones and the rate of the energetic in- 
stability to have an observable effect can be controlled 
by increasing the size of the thermal atom cloud 

A two-component ultra-cold ^^Rb vapor has also been 
trapped in a spin-dependent lattice using two counter- 
propagating laser beams with linear polarizations [3l|. 
The two species experience different <j+ and <j_ polarized 
optical lattices where the separation between the lattice 
potentials can be controlled by changing the angle be- 
tween the linear polarization vectors. 

The techniques developed for investigating dynam- 
ical and energetic instabilities in a single component 
case could be adapted to our proposed two-component 
BEC studies. A spin-dependent lattice potential may be 
used to control the value of the intraspecies interaction 
strength U12 by modifying the spatial overlap integral 
between the lattice site wavefunctions of the two species. 
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Moreover, the various intra- and inter-species scattering 
lengths and the two species Feshbach resonances [57] in 
two- component BEC systems make them a very rich area 
for experimental exploration. In addition to the two- 
component ^^Rb vapor [31], two-component BECs have 
been experimentall y r ealized in optical lattices using a 
'^^K-^'^Rb mixture [32^. 

Different superfluid velocities for the two species can be 
realized in such a system by moving the two lattice poten- 
tials at different speed. The disadvantage of this scheme 
is that it would make the intraspecies interaction strength 
Ui2 time-dependent. Displacing the harmonic trapping 
potentials of the two BECs different distances, in such a 
way that at the end the traps are perfectly overlapping, 
could be used to realize two BECs undergoing dipolar 
oscillations in phase with different amplitudes. Perhaps 
the easiest method to measure the reversed phase sepa- 
ration instability in a two-species BEC, as discussed in 
Section Bl B 31 and Appendix [Bf is to move the lattice po- 
tentials of the both species at the same speed and to use 
light-stimulated coherent Bragg diffraction [58] to change 
the velocity of one of the BECs. For sufficiently small ve- 
locities of the Bragg diffracted BEC, the atomic clouds 
of the two BECs overlap long enough for the dynamical 
instabilities to have an observable effect. 

Atomic ^'^Rb spin-1 gases have also been loaded to 
optical lattices [44]. In such a system the linear and 
quadratic Zeeman shifts could be modified, e.g., by us- 
ing off-resonant microwave field-induced level shifts [55]. 
This would allow the studies of the stability properties 
of different steady-state solutions presented here. 



V. CONCLUSIONS 

We studied the transport properties of two-component 
and spinor atomic BECs in optical lattices using the dis- 
crete nonlinear Schrodinger equations, obtained in the 
tight-binding approximation to the lattice system. The 
classical GP theory is valid in optical lattices at low tem- 
peratures if the effective ID nonlinearity is not too large, 
the atom number not too small, or the lattice potential 
not too deep [59|. In particular, we studied both the 
dynamical and energetic stability of homogenous Bloch 
wave solutions to the DNLSEs for the condensates by 
analyzing the linearized perturbations around the carrier 
wave. In the case of the dynamical instabilities this in- 
volved finding the eigenvalues (normal mode energies) of 
the corresponding Bogoliubov-de Gennes equations [the 
matrix aAi in Eq. ([5])] and in the case of energetic in- 
stabilities finding the eigenvalues of the second order per- 
turbations in the energy functional (the matrix Ai). Our 
steady-state Bloch wave ansatz allowed for magnetic Zee- 
man level shifts and even for two different velocities of 
the condensates in the two-component case. Our study 
discusses a large number of cases and points out how the 
spin degree of freedom can affect the stability properties. 

In the two-component case we analyzed and fully char- 



acterized the dynamical and energetic instabilities for the 
general case of the two BECs exhibiting arbitrary hop- 
ping amplitudes, interaction strengths, superfluid veloci- 
ties, and spatial overlap. Simple analytic expressions for 
the normal mode energies [Eq. (|6])] and the dynamical 
stability criteria were obtained in the important case of 
the two BECs having the same atom current (even when 
the velocities may differ). 

For the case that cos(/cia) and cos(/c2<^) exhibit equal 
sign for the two BEC carrier wavenumbers ki , k2 (see 
Sec. IIIB 2p . we found that the instability diagram con- 
tains contributions from: (1) the high velocity insta- 
bility, determined by Eq. ([12]), (which is analogous to 
that of a single-component BEC in a lattice), and (2) 
a weaker phase-separation instability [Eq. ([13])], which 
occurs when > ^11^22- However, as shown in Ap- 
pendix |B] and in Sec. Ill B 3\ an interesting case arises 
when one allows different condensate velocities of the 
two-components, so that cos(/cia) and cos(/c2<^) exhibit 
different sign (the effective masses of the two components 
exhibit different signs). Firstly, the high velocity instabil- 
ity conditions (which depends on the velocities, hopping 
amplitudes, atom numbers, and interaction strengths of 
the two BECs) indicate that the presence of the other 
condensate component can stabilize the superfluid flow 
of an otherwise unstable condensate (that exceeds the 
critical velocity of a single-component BEC). Secondly, 
the phase separation stability criteria can be reversed for 
particular sets of parameters and the entire dynamically 
stable regime exists for U12 > ^11^22; see Appendix [Bl 
and Figs.|3]^c-d). 

For the spin-1 BEC case, we also obtained analytic ex- 
pressions for the dynamical and energetic instabilities in 
several cases of interest. In the absence of the Zeeman 
level shifts the normal mode energies in the polar and 
ferromagnetic ground state manifolds are simplified and 
the two cases differ when Uq < and for relatively large 
spin-dependent interactions |/72| ~ |/7o|. In particular, 
the polar case tends to exhibit more regions of instability 
as Uq and U2 separately contribute [see Eq. ([T9|) ]. while 
in the ferromagnetic case, it is the sum Uq + U2 which is 
important; see Eq. ([25]) . This allows for dynamical sta- 
bility of ferromagnetic solutions for Uq < 0, and even for 
polar spin-dependent scattering lengths {U2 > 0). Also, 
we found that, unlike the polar case, the ferromagnetic 
solution will have some energetic instability for any fi- 
nite BEC velocity due to the existence of a pure kinetic 
energy eigenvalue ei± [Eq. ([26]) ]. 

In the presence of the linear and quadratic Zeeman 
level shifts we find a new set of steady-state Bloch wave 
solutions, describing the superfluid flow in spin-1 BECs. 
While the polar-like solution ([27|) and the ferromagnetic- 
like solution ([33]) form subsets of the corresponding so- 
lutions in the absence of the Zeeman splitting, this is 
not the case for the steady-state solutions (|3Q|) and (|34|) . 
The solution ([30]) exhibits a nonvanishing spin vector (F) 
pointing along the magnetic field and interpolates be- 
tween the polar and the ferromagnetic solutions. The 
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solution (|34|) only exists in the presence of sufficiently 
large Zeeman shifts and in this sense represents an en- 
tirely novel state. 

We analyzed the stability conditions for all spin-1 
Bloch wave states. For condensate solutions unique to 
the presence of Zeeman shift, Eqs. (|3Q|) and (|34j) . the 
stability diagrams were presented in Figs. [11] and [121 
respectively. For the parameters of ^^Na at low veloc- 
ities we found the solution Eq. (|3Q|) to be dynamically 
stable for negative quadratic shifts. Even for positive 
quadratic shift, a sufficiently large linear shift can sta- 
bilize it. Moreover, the solution (|3Q|) can be stable for 
ferromagnetic scattering coefficients {U2 < 0), even close 
to the polar state. For the parameters of ^^Rb the solu- 
tion Eq. (|34|) can be energetically and dynamically stable 
for positive quadratic Zeeman shifts and sufficiently large 
linear Zeeman shifts. 

The phenomena discussed in this paper should be ap- 
plicable to current and future experiments with multiple- 



component BECs in optical lattices, and we discussed 
some of the important considerations for the experimen- 
tal realization. We concentrated on the stability studies 
of moving Bloch wave solutions in the lattice. An in- 
teresting theoretical extension of this work is to consider 
inhomogenous condensate solutions, such as soliton-like 
structures in spinor BECs [47]. One could also investi- 
gate the effect of spin-dependent lattice potentials with 
spatially inhomogeneous profiles for the hopping ampli- 
tude, for instance dimerization, along the lattice [60]. 
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I 

APPENDIX A: BOGOLIUBOV-DE GENNES MATRICES 



In the two-component case, the Bogoliubov-de Gennes equation, upon substitution of the the ansatz (|4]) into ([2]), 
we get dU with: 



/ if 1+ + Ui 


-f/i 


U12 


-U12 \ 


-tJl 


Ki- + Ui 


-U12 


U12 


U12 


-U12 


K2+ + U2 


-U2 


V -U12 


U12 


-U2 


K2- + U2 J 



where Kj± = {AJj/h) sin^(ga/2) cos{kja) ± {2Jj/fi) sin(ga) smikjo)^ Uj = UjUj/h, and U12 = Ui2-sJnin2/Ti- 

In the general spinor case, we obtain the Ai for the spinor wavefunction and the Bogoliubov expansion (|4]), into 
the DNLSEs 



(^++)* m|+) (^0+)* (/oV)* (^-+)* (^-+)* 



M = 



.0 
Jo- 



^0+ 



5% 



^0+ 
.0- 

/o+ 
r-+ 



f°+ wo- M. 
{wo-r (/°+)* {w 



(-) 

-r 



Wo- 

/o- 

W-- 



(A2) 



M 



(+) 



(0) 



Mi") 

Jkl 

hi"- 

Wkl 



K± + U{1 + \c+f) + U{2\C+f + \Co\' - \C-f) + s+, 

K± + u{i + \Cof) + ui\c+f + \C-n 

K± + [7(1 + \C-f) + f7(-|C+P + \Cof + 2|C-n + S. 

{u - u)CkCi - uC, 

-UCl - 

-{u + u)CkCi, 
{u-u)CkCi. 



Here U = nUo/h, U = nU2/h, and 

K± = -2{J/n) cos{ka) 



{4:J/h) sm'^{qa/2) cos{ka) ± {2J/h) sm{qa) sm{ka) — ji . 



(A3) 
(A4) 
(A5) 
(A6) 
(A7) 
(A8) 
(A9) 
(AlO) 

(All) 
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APPENDIX B: DYNAMICAL STABILITY FOR 

kia > 7r/2, k2a < 7v/2 

In this section we analyze the dynamical stability of the 
two-component EEC system with equal atom currents 
Jisin(A:ia) = J2sin(/c2a), when cos(/cia) and cos (/c2a) 
exhibit different signs. Setting the atom currents to be 
equal allows us to obtain simple analytic expressions for 
the stability conditions. The different signs of cos (kia) 
and cos(A:2a) represent the situation where the veloci- 
ties of the two BECs are located on the opposite sides 
of the deflection point in the ideal, single-particle BEC 
excitation spectrum ([8|) (the effective masses of the two 
components exhibit different signs). Without loss of gen- 
erality we assume in the following that cos {kia) < and 
cos(A:2<^) > 0. The situation where when cos {kia) and 
cos {k2a) have the equal sign is covered in Sec IIIBl 

The analytic result for the normal mode energies is 
given by Eq. (|6]). Similarly to the case when cos {kia) 
and cos{k2a) have the same sign, the system is always 
dynamically unstable if Cjof 
same condition as in Eq, 



^2,q < we have the 



Dii + 1^22 < , 



(Bl) 



where Dij is defined in Eq. (pT]) . Since kia > 7r/2 and 
k2a < 7t/2 and Ji, J2 > 0, the inequality may even be 
satisfied for some values for which Uii^U22 > 0. In this 
case only one of the BECs reaches the (single-component) 
critical velocity ka = 7r/2, destabilizing the entire two- 
component BEC system. 

Next we assume ^i^g + ^2,g — ^ additional 
unstable regions of the parameter space. When cos {kia) 
and cos {k2a) have different signs, the expression inside 
the inner square root in Eq. (|6]) may become negative, 
resulting in a dynamical instability. In particular, this 
happens at least for some values of if 



where 



C^12 > a 



(Dn - £'22)= 



4nin2 Ji cos (/cia) J2 cos {k2a) 



(B2) 



(B3) 



where Dij is defined in Eq. (pT|) . 

Also the expression inside the outer square root may 
become negative. If 



niUii < 77.2^/22 min 



J2 cos{k2a) 
Ji cos{kia) 



Ji cos{kia) 
J2 cos(/c2a) 



(B4) 



this happens at least for some values of if U11U22 > 
Ui2' Combining this with Eq. (|B2p we find that the sys- 
tem is stable for the values of Uu that satisfy Eq. (|B4p . 
if 



U11U22 < < ^i . 
Similarly, for the values of Uu satisfying 



(B5) 



niUii > 712^/22 
niUii < n2U22 



Ji cos{kia) 



J2 cos{k2a) 
J2 cos{k2a) 



4Ji cos(A:ia) 



Ji cos{kia) 



(B6) 



we find that a dynamically stable system exists if 



(B7) 



where 



^2 =^11^22 +4Jicos(/cia)J2Cos(/c2a)/(nin2) 

+ 2[Ji cos{kia)U22 / ni + J2 cos{k2a)Uii/n2] • (B8) 

When Uu satisfies Eq. (jB6p . ^2 is always larger than 
U11U22 and for the dynamically stable region we have 
^11^22 < 6 < ^12 < ^^1- If ^2 > 65 no stable region 
exists. 



Note that the entire stable region in both Eqs. (jB5j) 
and (|B7p correspond to the values of the nonlinearities 
satisfying U11U22 < U12 that is normally associated with 
the dynamically unstable phase separation condition. 

The two-component system is therefore dynamically 
stable if Du + 1^22 > and U12 satisfies either Eq. (|B5p 
or Eq. dHll), for Uu defined by Eq. dHU) or Eq. re- 
spectively. Interestingly, we find a regime where the other 
condensate component can stabilize the superfluid flow 
of an otherwise unstable condensate. The inequalities 
(|B4[) or (|B6p can be satisfied for /7ii, U22 > when the 
component tjji exceeds the critical velocity of the single- 
component BEC, with kia > 7r/2, so that cc;^^ < in 
Eq. ([7j). The two-component BEC dynamics, neverthe- 
less, is stable if U12 satisfies either Eq. (|B5[) or Eq. (|B7j), 
respectively. 

Note also that we may have, e.g., Uu < 0, U22 > 0, 
but uf^q + cj2,g > (i.e., Du + 1^22 > 0). Such a two- 
component system can be dynamically stable since U12 > 

t/li/722. 
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